from numpy import *
import  matplotlib.pyplot as plt
def loadDataSet(url):
    with open(url) as fr:
        dataMat = [];labelMat = []
        for line in fr.readlines():
            lineArr = line.strip().split()
            dataMat = dataMat+[[float(lineArr[0]),float(lineArr[1])]]
            labelMat =labelMat+[float(lineArr[2])]
        return dataMat,labelMat

def selectJrand(i,m):
    j=i
    while(j==i):
        j = int(random.uniform(0,m))
    return j

def clipAlpha(aj,h,l):
    if aj>h:
        aj = h
    if l>aj:
        aj = l
    return aj

def smoSimple(dataMat,classLabels,c,toler,maxIter):
    dataMatrix =mat(dataMat)
    labelMat = mat(classLabels).T
    b= 0
    m,n = shape(dataMatrix)
    alphas = mat(zeros((m,1)))
    iters = 0
    while iters < maxIter :
        alphaPairsChanged = 0
        for i in range(m):
            fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T))+b
            Ei = fXi - labelMat[i]
            if (labelMat[i]*Ei<-toler and alphas[i]<c) or (labelMat[i]*Ei>toler and alphas[i]>0) :
                j = selectJrand(i,m)
                fXj = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T))+b
                Ej = fXj -labelMat[j]
                alphaIold = alphas[i].copy()
                alphaJold = alphas[j].copy() # numpy 复制要用copy()深度复制，不然会导致重新复制后，复制改变。这里可以理解为列表的函数操作
                if labelMat[i] !=labelMat[j] :
                    l = max(0,alphas[j]-alphas[i])
                    h = min(c,c+alphas[j]-alphas[i])
                else :
                    l = max(0,alphas[j]+alphas[i]-c)
                    h = min(c,alphas[j]+alphas[i])
                if l==h :
                    # print('l==h')
                    continue
                eta = 2* dataMatrix[i,:]*dataMatrix[j,:].T-dataMatrix[i,:]*dataMatrix[i,:].T-dataMatrix[j,:]*dataMatrix[j,:].T
                if eta>=0 :
                    print('eta>=0')
                    continue
                alphas[j] -= labelMat[j]*(Ei-Ej)/eta
                alphas[j] = clipAlpha(alphas[j],h,l)
                if abs(alphas[j]-alphaJold)<0.00001 :
                    # print('j not moving enough')
                    continue
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold-alphas[j])
                # alphai*yi+alphaj*yj = -k = alphai(new)*yi+alphaj(new)*yj
                b1 = b-Ei-labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T-labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].T
                b2 = b-Ej-labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T-labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j,:]*dataMatrix[j,:].T
                if 0<alphas[i]<c :
                    b = b1
                elif 0<alphas[j]<c :
                    b = b2
                else :
                    b = (b1+b2)/2
                alphaPairsChanged += 1
                # print('iters: %d i: %d, pairschanged %d' % (iters,i,alphaPairsChanged))
        if alphaPairsChanged == 0 :
            iters += 1
        else :
            iters = 0
        print('iteration number: %d' % (iters))
    return b,alphas

def dataPlot(dataArr,labelArr) :
    x1 = []
    y1 = []
    x2 = []
    y2 = []
    x = []
    y = []
    w = array((0.0, 0.0))
    for i in range(100):
        if alphas[i] > 0:
            x.append(dataArr[i][0])
            y.append(dataArr[i][1])
            w += array(alphas[i])[0][0] * labelArr[i] * array(dataArr[i])
    wx = arange(3, 6, 0.1)
    wy = -(array(b)[0][0] + w[0] * wx) / w[1]

    for i in range(len(labelArr)):
        if labelArr[i] > 0:
            x1.append(dataArr[i][0])
            y1.append(dataArr[i][1])
        else:
            x2.append(dataArr[i][0])
            y2.append(dataArr[i][1])
    plt.plot(wx, wy, color='r')
    plt.scatter(x, y, 32, color='r')
    plt.scatter(x1, y1, 12, color='g')
    plt.scatter(x2, y2, 12, color='b')
    plt.legend(['line','marker','-1','1'],loc=0)
    plt.show()

def kernelTrans(x,A,kTup) :
    m,n = shape(x)
    k = mat(zeros((m,1)))
    if kTup[0]=='lin' :
        k = x*A.T
    elif kTup[0] == 'rbf' :
        for j in range(m):
            deltaRow = x[j,:]-A
            k[j] = deltaRow*deltaRow.T
        k = exp(k/(-1*kTup[1]**2))
    else:
        raise NameError('Houston We Have a Problem That Kernel is not recognized')
    return k


class optStruct :
    def __init__(self,datain,classLabels,c,toler,kTup):
        self.x = mat(datain)
        self.labelMat = mat(classLabels).T
        self.c = c
        self.tol = toler
        self.m = shape(datain)[0]
        self.alphas = mat(zeros((self.m,1)))
        self.b = 0
        self.eCache = mat(zeros((self.m,2)))
        self.k = mat(zeros((self.m,self.m)))
        for i in range(self.m) :
            self.k[:,i] = kernelTrans(self.x,self.x[i,:],kTup)


def calcEK(os,k) :
    # fxk = multiply(os.alphas,os.labelMat).T*(os.x*os.x[k,:].T)+os.b
    fxk = multiply(os.alphas, os.labelMat).T * os.k[:, k] + os.b
    Ek = fxk-os.labelMat[k]
    return Ek

def selectJ(i,os,Ei) :
    maxk = -1
    maxDeltaE = 0
    Ej = 0
    os.eCache[i] = [1,Ei]
    validEcacheList = nonzero(array(os.eCache[:,0].A))[0]
    if len(validEcacheList)>1 :
        for k in validEcacheList:
            if k==i :
                continue
            Ek = calcEK(os,k)
            deltaE = abs(Ei-Ek)
            if deltaE>maxDeltaE :
                maxk = k
                maxDeltaE = deltaE
                Ej =Ek
        return maxk,Ej
    else:
        j = selectJrand(i,os.m)
        Ej = calcEK(os,j)
    return j, Ej

def updateEk(os,k) :
    Ek = calcEK(os,k)
    os.eCache[k] = [1,Ek]

def innerL(i,os) :
    Ei = calcEK(os,i)
    if (os.labelMat[i]*Ei<-os.tol and os.alphas[i]<os.c) or (os.labelMat[i]*Ei>os.tol and os.alphas[i]>0) :
        j,Ej = selectJ(i,os,Ei)
        alphaIold = os.alphas[i].copy()
        alphaJold = os.alphas[j].copy()
        if os.labelMat[i] != os.labelMat[j] :
            l = max(0,os.alphas[j]-os.alphas[i])
            h = min(os.c,os.c+os.alphas[j]-os.alphas[i])
        else:
            l = max(0,os.alphas[j]+os.alphas[i]-os.c)
            h = min(os.c,os.alphas[j]+os.alphas[i])
        if l==h:
            return 0
        # eta = 2*os.x[i,:]*os.x[j,:].T-os.x[i,:]*os.x[i,:].T-os.x[j,:]*os.x[j,:].T
        eta = 2*os.k[i,j]-os.k[i,i]-os.k[j,j]
        if eta >= 0:
            return 0
        os.alphas[j] -= os.labelMat[j]*(Ei-Ej)/eta
        os.alphas[j] = clipAlpha(os.alphas[j],h,l)
        if (abs(os.alphas[j]-alphaJold)<0.00001) :
            return 0
        os.alphas[i] += os.labelMat[i]*os.labelMat[j]*(alphaJold-os.alphas[j])
        updateEk(os,i)
        # b1 = os.b-Ei-os.labelMat[i]*(os.alphas[i]-alphaIold)*os.x[i,:]*os.x[i,:].T-os.labelMat[j]*(os.alphas[j]-alphaJold)*os.x[i,:]*os.x[j,:].T
        # b2 = os.b-Ej-os.labelMat[i]*(os.alphas[i]-alphaIold)*os.x[i,:]*os.x[j,:].T-os.labelMat[j]*(os.alphas[j]-alphaJold)*os.x[j,:]*os.x[j,:].T
        b1 = os.b - Ei - os.labelMat[i] * (os.alphas[i] - alphaIold) * os.k[i, i] - os.labelMat[j] * (os.alphas[j] - alphaJold) * os.k[i, j]
        b2 = os.b - Ej - os.labelMat[i] * (os.alphas[i] - alphaIold) * os.k[i, j] - os.labelMat[j] * (os.alphas[j] - alphaJold) * os.k[j, j]
        if 0<os.alphas[i]<os.c :
            os.b = b1
        elif 0<os.alphas[j]<os.c :
            os.b = b2
        else:
            os.b =(b1+b2)/2
        return 1
    else:
        return 0

def smop(datain,classLabels,c,toler,maxIter,kTup=('lin',0)) :
    os = optStruct(datain,classLabels,c,toler,kTup)
    iter = 0
    entireSet = True
    alphaPairsChanged = 0
    while iter<maxIter and (alphaPairsChanged>0 or entireSet) :
        alphaPairsChanged = 0
        if entireSet:
            for i in range(os.m) :
                alphaPairsChanged += innerL(i, os)
                # print('iters: %d i: %d, pairschanged %d' % (iters,i,alphaPairsChanged))
            iter += 1
        else:
            nonBoundIs = nonzero((os.alphas.A>0)*os.alphas.A<c)[0]
            for i in nonBoundIs:
                alphaPairsChanged += innerL(i, os)
                # print('iters: %d i: %d, pairschanged %d' % (iters,i,alphaPairsChanged))
            iter += 1
        if entireSet :
            entireSet = False
        elif alphaPairsChanged == 0:
            entireSet = True
        print('iteration number: %d' % (iter))
    return os.b,os.alphas

def testRbf(k1=1.3):
    fileurl = r'C:\Users\bgape\Desktop\computerS\machinelearninginaction\Ch06\\'
    dataArr,labelArr = loadDataSet(fileurl+'testSetRBF.txt')
    b,alphas = smop(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1)) #C=200 important
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    svInd=nonzero(alphas.A>0)[0]
    sVs=datMat[svInd] #get matrix of only support vectors
    labelSV = labelMat[svInd];
    print ("there are %d Support Vectors" % shape(sVs)[0])
    m,n = shape(datMat)
    errorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
        if sign(predict)!=sign(labelArr[i]): errorCount += 1
    print ("the training error rate is: %f" % (float(errorCount)/m))
    dataArr,labelArr = loadDataSet(fileurl+'testSetRBF2.txt')
    errorCount = 0
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    m,n = shape(datMat)
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
        if sign(predict)!=sign(labelArr[i]): errorCount += 1
    print ("the test error rate is: %f" % (float(errorCount)/m))

# fileurl = r'C:\Users\bgape\Desktop\computerS\machinelearninginaction\Ch06\testSet.txt'
# dataArr,labelArr = loadDataSet(fileurl)
# b,alphas = smop(dataArr,labelArr,1,0.001,40)
# dataPlot(dataArr,labelArr)
# print(b,alphas[alphas>0])

testRbf()


